First-order Chapman-Enskog velocity distribution function in a granular gas 
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A method is devised to measure tiie first-order Ciiapman-Enskog velocity distribution function 
associated witli tlie lieat flux in a dilute granular gas. The method is based on the application 
of a homogeneous, anisotropic velocity-dependent external force which produces heat flux in the 
absence of gradients. The form of the force is found under the condition that, in the linear response 
regime, the deviation of the velocity distribution function from that of the homogeneous cooling 
state obeys the same linear integral equation as the one derived from the conventional Chapman- 
Enskog expansion. The Direct Simulation Monte Carlo method is used to solve the corresponding 
Boltzmann equation and measure the dependence of the (modified) thermal conductivity on the 
coefficient of normal restitution a. Comparison with previous simulation data obtained from the 
Green-Kubo relations [Brey et al., J. Phys.: Condens. Matter 17, S2489 (2005)] shows an excellent 
agreement, both methods consistently showing that the first Sonine approximation dramatically 
overestimates the thermal conductivity for high inelasticity [a < 0.7). Since our method is tied to 
the Boltzmann equation, the results indicate that the failure of the first Sonine approximation is 
not due to velocity correlation effects absent in the Boltzmann framework. This is further confirmed 
by an analysis of the first-order Chapman-Enskog velocity distribution function and its three first 
Sonine coefficients obtained from the simulations. 

PACS numbers: 45.70.Mg, 47.57.Gc, 05.20.Dd, 51.10.+y 
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I. INTRODUCTION 

It is well established that granular fluids can be suc- 
cessfully described by fluid dynamics, nonequilibriuni 
statistical mechanics, and kinetic theory conveniently 
modified to account for the energy dissipation due to the 
inelasticity of collisions [1]. In particular, the Navier- 
Stokes (NS) transport coefficients characterize the de- 
parture from homogeneous situations due to small spa- 
tial gradients of the hydrodynamic fields (density, flow 
velocity, and granular temperature). In the case of a 
monodisperse granular fluid, the relevant transport coef- 
ficients are (i) the shear and the bulk viscosities (which 
relate the pressure tensor to the flow velocity gradients) 
and (ii) the thermal conductivity n and a coefficient /i not 
present for elastic collisions (which relate the heat flux to 
the temperature and density gradients, respectively). 

By using statistical-mechanical methods, formal ex- 
pressions for the transport coefficients in the form of 
Green-Kubo (GK) relations have been recently derived 
[2-7] . While these GK relations have a structure formally 
similar to their elastic counterparts, they are not simple 
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extensions of the latter and expose the subtleties associ- 
ated with the energy dissipation in collisions. In princi- 
ple, the GK formulae can be used to get the transport co- 
efficients by measuring the appropriate time correlation 
functions in the homogeneous cooling state (HCS) from 
computer simulations. In the low-density regime, and 
assuming the validity of the molecular chaos hypothe- 
sis, a kinetic theory description based on the Boltzmann 
equation is suitable [8, 9]. In that case, the GK rela- 
tions adopt a more explicit form, where the generator of 
dynamics is the linearized Boltzmann operator. These 
low-density GK expressions have been recently employed 
in computer simulations by means of the Direct Simula- 
tion Monte Carlo (DSMC) method to evaluate the depen- 
dence of the transport coefficients on inelasticity [10, 11]. 

In a kinetic theory description, an alternative route to 
determine the transport coefficients is provided by the 
Chapman-Enskog (CE) method [12]. Thus, the NS coef- 
ficients have been derived from the Boltzmann equation 
for monodisperse [13-15] and multicomponent [16] dilute 
gases, as well as from the Enskog equation for moderately 
dense gases [17, 18]. The equivalence between the GK 
and CE expressions derived from the Boltzmann equa- 
tion has been checked for several transport coefficients 
[3-5]. In the CE method, the transport coefficients are 
given in terms of the solutions of linear integral equations 
involving the linearized collision operator. Since the so- 
lutions are not exactly known, even in the elastic case, 
one usually resorts to the so-called Sonine approxima- 



tions. This allows one to get explicit expressions for the 
transport coefficients in terms of the parameters of the 
system. 

As in the elastic case, the first Sonine approxima- 
tion to the shear viscosity rj is seen to agree quite well, 
even for strong dissipation, with DSMC results obtained 
from three independent methods: (i) the time decay of a 
weak transverse shear wave in the HCS [19], (ii) the GK 
method [10, 11], and (iii) a modified simple shear flow 
[15, 20, 21]. Concerning the heat transport coefficients 
K and fi, a good agreement between the first Sonine pre- 
dictions and DSMC simulations of the corresponding GK 
formulae has been observed for values of the coefficient of 
normal restitution a > 0.7 [10, 11]. However, significant 
discrepancies appear for stronger dissipation {a < 0.7). 
While the range a > 0.7 widely covers those situations 
of practical or experimental interest, it is still important, 
from a fundamental point of view, to understand the ori- 
gin of those discrepancies for a < 0.7. This can shed 
light on the physical mechanisms playing a relevant role 
in the dynamics of granular flows at strong dissipation. 

Three scenarios are possible to explain the above dis- 
crepancies: (i) although the flrst Sonine approximation 
might accurately estimate the true coefficients n and 
fx described by the Boltzmann equation, the HCS ex- 
hibits relevant velocity correlations for strong inelastici- 
ties, even in the low-density limit, not accounted for by 
the inelastic Boltzmann equation; (ii) the discrepancies 
are an artifact of the DSMC method when applied to the 
evaluation of the time correlation functions needed in the 
GK representation; (iii) the disagreement is simply due 
to the limitations for high dissipation of the first Sonine 
approximation because of the deviation of the HCS ve- 
locity distribution from its Gaussian form. 

The aim of this paper is to contribute to the clarifi- 
cation of the above controversial issue. Specifically, we 
want to confirm or discard the possibility (iii) mentioned 
above. To that end, we will work within the framework of 
the Boltzmann equation and will use the DSMC method 
to compute one-time averages, in contrast to the GK 
route which involves two-time averages. Our method 
consists of perturbing the HCS by the application of a 
weak non-conservative (i.e., velocity-dependent) external 
force which produces a heat flux in the absence of in- 
homogeneities. This force, which mimics the effect of a 
thermal gradient, must be chosen under the condition 
that the perturbed velocity distribution function to first 
order coincides exactly with the one derived from the CE 
method at the NS order. This method is a non-trivial ex- 
tension to the granular case of the one devised time ago 
by Evans [22] and, independently, by Gillan and Dixon 
[23]. However, while in the elastic case the force is par- 
allel to the heat flux, this is not so in the inelastic case, 
as a consequence of the non-Gaussian form of the HCS. 
For elastic systems, the method has been successfully ap- 
plied to get the thermal conductivity of Lennard-Jones 
[24] and hard-sphere [25] fluids; the corresponding non- 
linear problem has also been studied [26, 27]. Apart from 



measuring the transport coefficients, the method allows 
one to get the velocity distribution function to first order 
in the external field and this is one of the main objectives 
of this paper. 

As will be seen later on, in order to determine the 
thermal conductivity k and its associated velocity distri- 
bution, the application of an external force is not enough 
since an additional stochastic term is needed, which com- 
plicates the implementation of the simulation method. 
The situation is even worse in the case of ^ since it is 
coupled to K. Nevertheless, the integral equation associ- 
ated with the coefficient k' = k — nfi/2T can actually be 
simulated by the action of an external force only. There- 
fore, in this paper we focus on this modified thermal 
conductivity coefficient k' and its corresponding veloc- 
ity distribution. The simulation results obtained from 
the present method for k' agree well with those obtained 
from the alternative GK method [11]. Moreover, the ve- 
locity distribution function for high inelasticity differs ap- 
preciably from the Sonine expansion truncated after the 
first term. These results strongly support the scenario 
(iii) mentioned before. In fact, a modified first Sonine 
approximation, where the Gaussian weight appearing in 
the Sonine expansion is replaced by the HCS distribu- 
tion, compares fairly well with the simulation data for 
the whole range of inelasticities [28] . 

The organization of the paper is as follows. The Boltz- 
mann equation for inelastic hard spheres and its solution 
provided by the CE method is presented in Sec. H. The 
linear integral equations for the NS distributions associ- 
ated with the heat flux are worked out and the modifled 
thermal conductivity coefficient k' is introduced. In Sec. 
Ill we propose the method to get the first-order (NS) dis- 
tributions by the application of velocity-dependent ex- 
ternal forces in spatially uniform states. The numeri- 
cal results obtained from the DSMC method for k' and 
its corresponding velocity distribution function are pre- 
sented and discussed in Sec. IV. The paper is closed with 
the conclusions in Sec. V. 



II. INELASTIC BOLTZMANN EQUATION AND 
CHAPMAN-ENSKOG EXPANSION 

For the sake of completeness and to fix the notation, 
in this Section we summarize the main known results 
obtained by applying the CE method to the Boltzmann 
equation. 

Let us consider a granular gas composed by smooth 
inelastic disks {d ~ 2) or spheres {d = 3) of mass m and 
diameter a. The inelasticity of collisions among all pairs 
is characterized by a constant coefficient of normal resti- 
tution a < 1. In the low-density regime, the evolution of 
the one-particle velocity distribution function /(r, v, t) is 
given by the Boltzmann kinetic equation [8, 9]. In the ab- 
sence of external forces, this equation reads, in standard 
notation, 



(9t + vV)/(r,v;t) = J[v|/,/], 



(2.1) 



where J[v|/, /] is the Boltzmann cohision operator [8, 9]. 
Cohisions conserve mass and monicntuni, but energy is 
dissipated (except, of course, if a = 1): 



/dvJ V Ij[v|/,/] = J 



m 3 



Here 



n = / dv/(v) 



(2.2) 



(2.3) 



is the number density, V = v — u is the pecuhar velocity, 
where 



u-i /dv/(v) 
n J 



is the flow velocity, 



T 



m 
dn 



d^fV'^fiw) 



(2.4) 



(2.5) 



is the granular temperature, and Q is the cooling rate. 
Similarly, the fluxes can be obtained as moments of the 
velocity distribution function. In particular, the irre- 
versible part of the pressure tensor is 



n, 



m I d-v { ViVj 



^%, 



/(v) 



and the heat flux is 



m 



j dwV^YJ{w). 



(2.6) 



(2.7) 



The CE method assumes the existence of a normal so- 
lution in which all the space and time dependence of the 
distribution function appears through a functional de- 
pendence on the hydrodynamic fields. For small spatial 
variations, this functional dependence can be made local 
in space and time through an expansion in gradients of 
the fields: / = Z^"-* -I- /^^^ -!-•••. The local reference state 
/^°' is constrained to have the same first few moments 
(2.3)-(2.5) as the exact distribution /. The kinetic equa- 
tion for /(o) is [13] 

ic("^^-(v/("))=J[V|/(°),/(°)], (2.8) 

where C'-^^ is defined by setting / -^ /*^°) in Eq. (2.2). 
Equation (2.8) can be recognized as the one satisfied by 
the HCS [8, 29], parameterized by the local values of n, 
u, and T. The exact solution to Eq. (2.8) has not been 
found, although some of its properties are known [29, 30]. 
Since /'"'(V) is an isotropic function, its Sonine ex- 
pansion is 



/(°)(V)-/m(V) 



l + Y^a,L[-^\c^) 



fc=2 



(2.9) 



where 



/m(V)=™o-''^ 



-d-d/2-c^ 



is the Maxwellian distribution, 

V 



vo{T) 



(2.10) 



(2.11) 



is the velocity relative to the thermal speed vo{T) = 
y^2T/m and Lj^ (x) are the generalized Laguerre poly- 
nomials whose orthogonality relation is [31] 



' dxxPe--LiP\x)L^^\x) = £i^±i±Z)4,. (2.12) 



The coefficients a^ are related to the velocity moments, 
measuring the deviation of /("' from the Gaussian. In 
particular, 02 is the fourth cumulant (or kurtosis) of the 
distribution function f^^': 



j dVV^f^\Y) =d{d + 2 



^(1 + 02). 



(2.13) 



By inserting Eq. (2.9) into Eq. (2.2) and neglecting ak 
with fc > 3 and nonlinear terms in 02, one gets 



^(0) 






where 



TT^ 



fo 



d + 2 r(rf/2) "'^ \mj 



1/2 



(2.14) 



(2.15) 



is an effective collision frequency. An excellent estimate 
of 02 is [32, 33] 



0-2 



16(l-a)(l-2Q;2) 



25 + 24d - a(57 - 8d) - 2(1 - a)a^ ' 



(2.16) 



To first order, the application of the CE method yields 
the linear integral equation [13] 



af^+/:)/(i)(V) =A(V)-Vlnr + B(V)-Vlnn, 



(2.17) 
where we have particularized to the case ViUj = 0. In 

Eq. (2.17), 9( is an operator acting on any function of 
temperature as 

di"^XiT) = {0i"h) ^ = -C("^r^, (2.18) 

C is the linearized Boltzmann collision operator defined 
as 



CX{V) = -J[V\,r>,X]-J[V\XJ^^>] 



(0)1 



(2.19) 



and 

A(V) 



2 dV 



V/(o)(V) 



T d 
m 9V 



/(°)(V), (2.20) 



B(V) ^ -V/(")(V) - - A/(o)(v). (2.21) 

m o\ 

The structure of the solution to Eq. (2.17) is 

/(^) (V) = ^(V) • V In T + B(V) • V In n. (2.22) 

Taking into account Eq. (2.7), the heat flux to first order 
(NS order) is 



.(1) - 



-kS/T — fj,Vn, 



where 



-^|dVS(V).AV) 



is the thermal conductivity coefficient and 



(2.23) 
(2.24) 

(2.25) 



is a new coefficient absent for normal gases. In Eqs. (2.24) 
and (2.25) we have introduced the function 



S(V) = 



-V'- 



d + 2^ 



T V. 



(2.26) 



By dimensional analysis, »4.(V) — nv^ A^*(c), where 
A ~ 1/na is the mean free path and ^*(c) is a dimen- 
sionless function of the reduced velocity c defined in Eq. 
(2.11). A similar relation holds for iB(V). Consequently, 



9(")/-(i)(v) = i^wA 



v/(i)(v) 



-C^°)A(V) • (vinn+-VlnT 

(2.27) 

Equating the coefficients of V In T and V In n in Eq. 
(2.17), one obtains the following pair of linear integral 
equations: 

C + ic<"^ ^ • V - IC^"A A{V) = A{V), (2.28) 



C + -C^"^S- ■ V^ ^(V) - C^"U(V) = B{V). (2.29) 
2 av / 



While Eq. (2.28) is a closed equation for the unknown 
fimction A, Eq. (2.29) is coupled to Eq. (2.28), so that 
one needs first to know A. to determine J3. However, the 
combination A! = A— ^B verifies the closed integral 
equation 



where 



c + Ic^"^^-v]a'{y) = a'{y), 



A'{V) ^ A(V) - -B(V) 



(2.30) 



d_ 



G'(V)/(o)(V) 



(2.31) 



and we have called 



T 



G'.^)^^^^-^r 



(2.32) 



Thus, the function A'(V) can be expressed as the diver- 
gence of a tensor. Let us see that the same applies to the 
function A{V). Note first the identity 



d 



av 



= -Vf(")^ 



Yf'\Y) =-v/^"^(v) 



d 

dV 



vv/(")(v) 

(2.33) 

Further, taking into account that f^^'(V) is an isotropic 
function, we can write 



V 



av' 



V/(")(V)] = (d-2)V/(")(V) + ^ [v\f^°\V) 



(2.34) 

Equating the right-hand sides of Eqs. (2.33) and (2.34), 
one has 



T/,/(°)(V) 



1 



d 



d - 1 dV, 



(T/.y,-y%)/(")(v) 



Insertion of this into Eq. (2.20) yields 



d 



A(V) = -— .[G(V)/(")(V) 



(2.35) 
(2.36) 



where the tensor G(V) is 



G.,(V) = - 



V' 



T 



2{d-l) m 



5ii - 



d-2 



j^^V.Vr (2.37) 



The problem of determining the functions ^(V) and 
5(V) from Eqs. (2.28) and (2.29) is fully equivalent to 
that of determining the functions ^(V) and ^'(V) from 
Eqs. (2.28) and (2.30), respectively. The latter wiU be the 
point of view adopted in the remainder of this Section as 
well as in Sec. III. In particular, the coefficient ^ defined 
by Eq. (2.25) can be expressed as 



IT 

— (k-k ), 
n 



where 



.±,JdVSiV)-A\V) 



(2.38) 



(2.39) 



is a modified thermal conductivity coefficient. In terms 
of it, Eq. (2.23) can be rewritten as 



.(1) 



-k'VT - ^r-l/2y friT^/2 



(2.40) 



Therefore, one has q'^-* = — kVT in those points where 
the density gradient vanishes, while q'^' — —k'WT in 
those points where the gradient of the cooling rate (which 
is proportional to nT^'"^) vanishes. In this context, it is 
worth mentioning that one of the hydrodynamic modes of 
a dilute granular gas corresponds to an excitation where 



(■(0) — const at zero flow velocity [34]. In the clastic case, 
/(o) is the Gaussian, C^°^ = 0, and A(V) = A'(V), so 
that K. ~ k' and ^ — Q. 

Formal expressions for k and k' can be derived by mul- 
tiplying both sides of Eqs. (2.28) and (2.30) by S(V) and 
integrating over velocity. The results are 



(2.41) 



d + 2nT 


l + 2a2 


2 m 


i^. - 2C(") ' 


d + 2nT 


l + |a2 



m Vk 



- f C^"^ ' 



(2.42) 



where 



j d\ S(\) ■ CA{V) 

/dvs(v).^(v) ■ 



_ /dVS(V)-£A^(V) 

"' " /dvs(v)-A'(v) ■ 

(2.43) 
The above expressions are formally exact, but the depen- 
dence of J^K, Vk.1 , C}-'^\ and 02 on a is unknown. While the 
two latter quantities require the knowledge of the HCS 
distribution f^ \ the collision frequencies v^ and v^i are 
given in terms of the solutions of the two linear integral 
equations (2.28) and (2.30). As said before, C^°^ and 02 
can be well approximated by Eqs. (2.14) and (2.16), re- 
spectively. 

The symmetry properties of the vectorial functions 
^(V) and ^'(V) suggest the following Sonine expansion 
representations : 



A{Y) 
A{Y) 



fc=i 



fMmi:\y\L'r\<i 



(2.44) 



Making use of the orthogonality condition (2.12), the co- 
efficients bk and b'f, can be expressed as moments of A. 
and A , respectively. In particular, the first coefficients 
&i and b'l are directly related to the thermal conductivi- 
ties as 



d + 2 



nXvo 



bi 



(2.45) 



Unfortunately, when the expansions (2.44) are inserted 
into Eqs. (2.28) and (2.30), one gets an infinite hierarchy 
of equations for the coefficients bk and b'j^, so that k and k' 
cannot be obtained exactly. For practical purposes, it is 
usual to truncate the Sonine expansions at a given order 
k and solve an approximate set of k equations. The case 
fc = 1 yields the so-called first Sonine approximation. In 
this case, the collision frequencies defined in Eq. (2.43) 
are [14] 



fo- 



1-t-a 



d — 1 3 , , „, ,, , 

-^ + -{d + 8){l-a) 



4 + 5d-3 4-d a 

H ^ —a2 

512 



(2.46) 



Insertion of Eqs. (2.14), (2.16), and (2.46) into Eqs. (2.41) 
and (2.42) gives the transport coefficients k and k' in the 



first Sonine approximation. The coefficient /i is then ob- 
tained from Eq. (2.38). In the elastic case (a = 1), both 
K and k' reduce to the thermal conductivity coefficient 
for elastic hard spheres, which in the first Sonine approx- 
imation is given by 



KO = 



d{d + 2) nT 
2{d ~ 1) ravQ 



(2.47) 



III. HOMOGENEOUS STEADY HEAT FLOW 

DRIVEN BY A VELOCITY-DEPENDENT 

EXTERNAL FORCE 

Apparently, the most direct way of measuring the NS 
velocity distribution function (2.22) as well as its asso- 
ciated transport coefficients k. and k! would imply the 
introduction of weak temperature and density gradients. 
However, this gives rise to several technical problems. 
First, the coupling between inelasticity and spatial gra- 
dients make it difficult to extract the real NS contribu- 
tions. Moreover, one must be able to identify the bulk 
region, where the boundary effects are negligible. In ad- 
dition, the measured quantities are local and hence the 
statistical errors become significant. Finally, it would 
be difficult to disentangle the contributions to the heat 
flux coming from the temperature and density gradients. 
This encourages the search for alternative methods that 
avoid these difficulties. 

As said in the Introduction, we want to study a ho- 
mogeneous nonequilibrium steady state generated by the 
action of an anisotropic velocity-dependent external force 
of strength e, which induces a heat flux in the absence 
of any thermal and density gradients. The form of the 
force must be chosen under the condition that, in the 
limit e ^ 0, the deviation of the velocity distribution 
function from that of the HCS is the same as the one 
produced by real thermal and density gradients. In the 
latter situation, the deviation from the HCS is measured 
by the NS functions -4.(V) and ^'(V), as discussed in 
the preceding Section. Therefore, we have to deal with 
two separate homogeneous steady Boltzmann equations, 
each one reducing in the linear order in e to Eqs. (2.28) 
and (2.30), respectively. 



A. Navier— Stokes function ^'(V) 

Let us start with the function .4.'(V) associated with 
the modified thermal conductivity k' . Taking into ac- 
count the structure of Eq. (2.30), together with Eq. 
(2.31), we propose the following Boltzmann equation: 






[F'(V)/(V)]-J[V|/,/] 



(3.1) 



where F'(V) is an external force (except for a factor m). 
More specifically, F' = F'^ -I- Fth is decomposed into a 



heat-flux force 



F',(V) = G'(V).e 



and a thermostat force 



Fth(V) = |V + /3. 



(3.2) 



(3.3) 



In Eq. (3.2), the tensor G'(V) is given by Eq. (2.32) and 
e is a constant vector which mimics the effect of V In T 
(at VnT^'^ = 0) and whose magnitude e measures the 
strength of the force. The parameters £, and (3 in the 
thermostat term are introduced to keep the total kinetic 
energy and momentum constant, respectively. 

Multiplying both sides of Eq. (3.1) by V, integrating 
over velocity, and imposing a vanishing total momentum 
we get the following expression for the parameter (3: 



(3^ 



2mn 



-n. 



(3.4) 



where Hij is defined by Eq. (2.6). Analogously, the con- 
dition of constant total kinetic energy yields 



^ = C + 



dnT 



qe. 



(3.5) 



where the cooling rate ( and the heat flux q are defined 
by Eqs. (2.2) and (2.7), respectively. Since the direction 
of e can be chosen arbitrarily without loss of generality, 
in what follows we take e = ex. By symmetry reasons, 
the only non-zero elements of Hij are Hxx and Ilyy = 
Ilzz = ■■■ = ^dd- Analogously, q = q^^ and /3 = P^^. In 
addition, in the limit e ^ 0, one has the leading behaviors 
qx ^ £ and Hij ^ e^, so that f3x ~ e^ and C ~ C ^ £^- 

In principle, Eq. (3.1) is a highly nonlinear Boltzmann 
equation since not only the collision term is quadratic 
in /, but also the coefficients ^ and /3 are functionals 
of / through (, Hij, and q. The control parameter in 
Eq. (3.1) is the field strength e. Although the problem of 
solving Eq. (3.1) for finite e is interesting by itself [26, 27], 
here we focus on the regime of small e. In that case, the 
solution to Eq. (3.1) can be expanded in powers of e as 



/(V) = /(°)(V) + /(i)(V) + 0(62). 



(3.6) 



Note that the expansion (3.6) does not need to be conver- 
gent but only asymptotic, similarly to what happens in 
the CE expansion. Setting e = in Eq. (3.1) and taking 
into account that ^ = C^") -I- O(e^), it is straightforward 
to see that /'•^•'(V) verifies Eq. (2.8), i.e., the Boltzmann 
equation for the HCS. Next, to first order in e, one has 



d 



F',(v)/(")(V) + ic*°^v/(^)(v) 



Therefore, 



/(i)(V)=A'(V).e, 



-£/(i)(V). 
(3.7) 

(3.8) 



where -4.'(V) obeys the linear integral equation (2.30) 
and use is made of Eq. (2.31). This proves that the depar- 
ture from the HCS distribution produced by the external 
force F'(V) coincides to first order in e with the one pro- 
duced by a real thermal gradient. As a consequence, the 
modified thermal conductivity k' can be obtained from 
the solution to Eq. (3.1) through the linear response re- 
lation 



/ 1 ,. fe 

K = --lim — 

1 e^O e 



(3.9) 



Analogously, the coefficients b'f. defined in Eq. (2.44) can 
be evaluated in terms of averages of velocity polynomials: 



b'. 



where 



2ni + d/2)k\ {cxLf'\^)) 

r(fc-fi-hd/2)e.™o 7* ' ^^-^"^ 

{X{V)) = - [dVXiV)f{V) (3.11) 

n J 



and e* = Ae is the reduced field strength, which plays the 
role of a Knudsen number. Since /'"'(V) is an isotropic 
function, only /'^•'(V) contributes to the averages of Eq. 
(3.10). Furthermore, it is proven in the Appendix that 
the coefficients feJc can alternatively be evaluated as 



K- 



2r(3/2)fc! 



(1/2) 



r(fc 



, , lim ■ 

3/2) e*^0 
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(3.12) 



While the averages in Eq. (3.10) require the knowledge of 
the full distribution /'^'(V), the averages in Eq. (3.12) 
only require the knowledge of the marginal distribution 
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g''>iVx)= dV^r^'iV 



Wi 



(3.13) 



where V^ = V — VxX. The equivalence between Eqs. 
(3.10) and (3.12) implies that all the information con- 
tained in the full distribution /^^'(V) is encapsulated in 
the marginal distribution g^^'{Vx). 

For sufficiently small values of e, Eq. (3.6) can be trun- 
cated, so that / — /(°) 4- f'-^\ At the level of the corre- 
sponding marginal distributions. 



g{Vx)^9^"Hv.)+9^'Hv.). 



(3.14) 



By symmetry, g^'^'^iVx) is an even function of Vx, while 
g^^'^iVx) is an odd function. Consequently, 



9^'HVx)^T;[9iV.)-g{-Vx)]. 



(3.15) 



This relation is useful for extracting the first-order dis- 
tribution g^^' (Vx) from the complete distribution g{Vx), 
provided that e is small enough to neglect nonlinear 
terms. 



B. Navier— Stokes function A.(V) 

Comparison between the integral equations (2.28) and 
(2.30) shows two differences. First, the inhomogeneous 
terms are different, i.e., A(V) 7^ A'(V). This means 
that the role of the force (3.2) is now played by the force 



F,(V) = G(V) • e. 



(3.16) 



where G(V) is given by Eq. (2.37). The most significant 
difference is the presence of the extra term — ^(^^^^ Ji(V) 
on the left-hand side of Eq. (2.28). This term cannot 
be accounted for by an external force. Therefore, the 
appropriate homogeneous steady Boltzmann equation in 
this case is 



^ • [F(V)/(V)] 



J[V\f,.f] + lc{.f-.f^"^), (3.17) 



with F = Fe + Fth, where the heat-flux force is given by 
Eq. (3.16) and the thermostat force is again (3.3). The 
condition of zero total momentum yields 



A 



d-2 



2{d — l)mn 



n, 



(3.18) 



while, by the condition of constant energy, ^ is still given 
by Eq. (3.5). Again, we can choose e — ex. As before, 
in the limit e ^ 0, one has qx ~ e, n^j ~ e^, f3x ^ e'^, 
and S, — C ^ ^ ■ Inserting the expansion (3.6) into Eq. 
(3.17), and following the same steps as in the preceding 
subsection, one sees that /(°) is the HCS distribution, 
Eq. (2.8), and 



/W(V)==A(V).e, 



(3.19) 



where -4.(V) obeys the linear integral equation (2.28). 
Thus, the thermal conductivity k can be obtained from 
the solution to Eq. (3.17) through a Hnear response re- 
lation similar to Eq. (3.9). Analogously, the Sonine co- 
efficients &fc can be evaluated from expressions similar to 
Eqs. (3.10) and (3.12). 

The main peculiarity of the Boltzmann equation (3.17) 
lies in the presence of the anomalous term ^C (/ — /''°-') 
on the right-hand side. This term represents a stochastic 
donation-annihilation process. According to it, a parti- 
cle of a given velocity V has a probability per unit time 
^C of being "cloned" and a probability per unit time 
1(^/(0) (V)//(V) of being "annihilated." Since / and 
f^^' are normalized to the same number density, flow ve- 
locity, and temperature, the donation-annihilation pro- 
cess does not change the total number of particles, mo- 
mentum, and energy. However, this stochastic process 
tends to increase the departure of the distribution func- 
tion from that of the HCS. It is worth remarking that, 
in the homogeneous steady Boltzmann equation associ- 
ated with the shear viscosity [21], a similar stochastic 
process is also present, but with the opposite sign, so 
that the stochastic term tends to decrease the departure 
from of the HCS. The presence or absence of this type 



of stochastic term has a counterpart in the GK relations. 
While the time correlation function is multiplied by an 
exponentially growing factor in the case of k, there is an 
exponentially decaying factor for 77, and there is no factor 
in the case of k' [5] . 

The presence of the non-standard term ^( (/ — /^"') 
makes it difficult, but not impossible, to implement the 
Boltzmann equation (3.17) by the DSMC method. Thus, 
for the sake of simplicity, henceforth we will focus on the 
Boltzmann equation (3.1) to determine the NS distribu- 
tion function ^(V) and its associated transport coeffi- 
cient k'. 

Before closing this Section, it is worth mentioning that 
the heat-flux forces (3.2) and (3.16) are not straightfor- 
ward extensions of the one proposed by Evans [22] and 
Gillan and Dixon [23] for normal fluids. The latter is 
Ff{V) = G'='(V) • e, where 



G^(V) 



dT 
2m 



V 



(3.20) 



As a consequence, F°'(V)||e, while FjCV) and F^(V) are 
not paralld to e. For a < 1, ^v • [G°'(V)/(")(V)] ^ 
dv ■ [G(V)/(")(V)] ^ dv ■ [G'(V)/(o)(V)]. However, in 
the elastic case (a = 1), /'•^•'(V) = /m(V), so that 



d 



[G^'(V)/m(V)] 



o> 

dV 
d 



[G(V)/m(V)] 
[G'(V)/m(V)] 



= r-iS(V)/M(V), (3.21) 
where S(V) is given by Eq. (2.26). 



IV. RESULTS 

In this Section we present results obtained by numeri- 
cally solving the Boltzmann equation (3.1) in the three- 
dimensional case {d — 3) by means of the DSMC method 
[35]. It has been rigorously proven that the DSMC 
method produces a solution to the Boltzmann equation in 
the limit of vanishing discretization and stochastic errors 
[36]. Details of the method can be found elsewhere [35], 
so that here we only point out some specific aspects in 
our simulations. The system consists of N simulated par- 
ticles, whose velocities are updated from time t to time 
t + 5t due to (i) the collisions and (ii) the action of the 
external force F'(V). The collisional stage (i) proceeds 
as usual [35], except that the collisions are inelastic. The 
streaming stage (ii) proceeds in two steps: 

1. Update the velocity of every particle i = 1, . . . ,N 
as follows: 

V^(t + St) = Ve{t) - iv,(t)V,(t) • eSt. (4.1) 

This only accounts for the velocity-dependent part 
of the heat-flux term F^(V), as can be seen from 
Eq. (2.32). 



2. Compute the mean velocity u'{t + St) 
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N ^J2e=i'^ei^ + ^^) ^^^ ^^^ kinetic energy per 
particle K'{t + 5t) = N-\m/2) X;^i [V^(i + St) - 
u'{t + St)]'^. In general, u'(i + St) ^ because 
of the step (4.1) and K'{t + 5t) ^ Kq, where Kq 
is the prescribed initial kinetic energy per particle, 
because of the collisional stage and also because of 
the step (4.1). Next, update again the velocity of 
every particle as 



^i{t + St) 



Kn 



K'{t + St) 



[Y',{t + St)-u'{t + St)]^ 



(4.2) 
This change guarantees that the total momentum 
is restored to zero and the kinetic energy to its ini- 
tial value, before proceeding again to the collisional 
stage. Equation (4.2) incorporates the effect of Fth, 
as well as of the velocity-independent part of F^ . 

Once a steady state is reached, the most relevant quan- 
tities can be computed. This includes the heat flux 

qx, the averages {cxL). (c^)) and {cxL\. (c^)) with 
k — 1,2,3, and the marginal velocity distribution func- 
tion g{Vx)- Note that k' and h'^ are essentially the same 
quantity, as shown by Eq. (2.45). By assuming that the 
value of the reduced strength e* is small enough to probe 
the linear regime, the modified thermal conductivity k! 
and the Sonine coefficients b'l, b'2, and 63 are obtained 
from Eqs. (3.9), (3.10), and (3.12). In addition, the 
marginal NS distribution function g^^'{Vx) is obtained 
from g{Vx) by applying Eq. (3.15). 

The number of simulated particles is A^ = 2 x 10^ and 
the time step is St = 0.003r, where r — A/wo is the 
mean free time and A — l/\/27m(T^. To improve statis- 
tics, the results are averaged over 200 independent re- 
alizations. The data are further averaged over time in 
the steady regime. The range of inelasticities analyzed 
is 0.3 < a < 1 and in all the cases the initial state is a 
Gaussian distribution. 

We first analyze the evolution toward the steady state. 
In the transient regime one can measure effective time- 
dependent values of k' and b'^^. Figure 1 shows the ratio 
K'{t)/Ko, where kq is defined by Eq. (2.47), versus t/r 
for a = 0.3 and two different values of the reduced field 
strength, namely e* = 0.05 and e* = 0.025. Even in this 
highly inelastic case, we observe that a steady-state value 
of k' is reached after about ten collisions per particle. 
Also, the figure clearly indicates that the transient and 
the stationary results are very weakly dependent of e*, 
i.e., the heat fiux Qx is practically linear in e*. We have 
observed a similar behavior for other values of a. Thus, 
we take the value e* = 0.025 in the remainder of the 
figures. An even smaller value of e* would decrease the 
signal-to-noise ratio without affecting the results. 

An additional test of the linear regime is provided by 
Fig. 2, which shows the time evolution of the Sonine coef- 
ficients b'l, b'2, and 63 for a = 1 and a = 0.3. These coeffi- 
cients have been evaluated from Eq. (3.10) and also from 
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FIG. 1: (Color online) Plot of the time-dependent (modified) 
thermal conductivity versus time for a — 0.3 and e* — 0.05 
(dashed line) and e* — 0.025 (solid line). The horizontal 
dotted line indicates the steady-state value averaged over time 
in the case e* — 0.025. 



Eq. (3.12). As proven in the Appendix, both methods 
must lead to identical results, provided that the system 
is in the linear regime. Figure 2 shows that the simulation 

data obtained from [cxL], (c^)) are practically indistin- 
guishable from those obtained from {cxL\, (c^)), thus 
confirming that the value e* = 0.025 is small enough to 
assume that the system is actually in the linear regime. 
Given that the fluctuations are smaller when the Sonine 
coefficients are computed from {cxL\. (c^)) than when 

they are computed from {cxL\. (c^)), henceforth we 
choose the former method. We also see in Fig. 2 that 
the characteristic relaxation time {t/r ~ 10) is practi- 
cally independent of the degree of inelasticity as well as 
of the order of the Sonine coefficient considered. 

Now we focus on the dependence of the steady-state 
quantities on dissipation. Figure 3 shows k'/kq as a func- 
tion of the coefficient of restitution a. In addition to our 
simulation data, we have included the values of k' ob- 
tained from the simulation data of k and /i reported by 
Brey et al. [11], as well as the first Sonine approximation 
given by Eqs. (2.42) and (2.46). It must be emphasized 
that the methods used here and in Ref. [11] are quite dif- 
ferent. In the latter method, two-time correlation func- 
tions of certain quantities are evaluated in the HGS and 
subsequently those correlation functions are integrated in 
time to get the transport coefficients. On the other hand, 
in our method the HCS is slightly perturbed by a homo- 
geneous external forcing and the transport coefficient is 
determined from a standard one-time average (namely, 
the heat flux), by assuming linear response. While the 
first Sonine prediction does a good job for a > 0.7, it 
dramatically overestimates the thermal conductivity for 
higher inelasticity. It is important to remark that the ex- 
cellent agreement found between both simulation meth- 
ods strongly supports the conclusion that the discrepan- 
cies between the first Sonine approximation and simula- 
tions are not due to the presence of velocity correlations 
beyond the Boltzmann description for strong dissipation. 



z,.u 


' 1 ' 1 


' 1 


1 ' 


1 ' 


1.5 




^3 

\ 






- 






a=i ; 


-^■^0.5 




- 


0.0 


1 1 1 1 


1 1 


1 1 


1 1 


2 
1 


' 1 ' 1 


' 1 


1 ' 


1 ' 






jtfciejsW^Oi. 


t'^^t^i^jtO'-i^ 


-1 


a=0.3 


■ 




1 1 1 1 


1 1 


1 1 


1 1 



1 St Sonine approxim. 

A GK 

• This work 



5 10 15 20 25 30 

tk 

FIG. 2: (Color online) Plot of the time-dependent Sonine co- 
efficients &J. (k = 1, 2, 3) versus time for a = 1 (top panel) and 
a — 0.3 (bottom panel). The solid and dahed lines correspond 
to data obtained from {cxL^, (c^)) and from {cxLl, (c^)), 
respectively. Note that they are distinguishable only in the 
case of &3 for a — 0.3. The horizontal dotted lines indicate 
the steady-state values. 



Instead, those discrepancies are essentially due to the in- 
accuracy of the first Sonine approximation for a < 0.7. In 
the elastic limit, the numerical results show that the first 
Sonine approximation slightly underestimates the ther- 
mal conductivity. This is known to be corrected if higher 
orders in the Sonine polynomial expansion are taken into 
account [12], which yields k/kq = 1.025218 [37]. This 
value is indicated by an arrow in Fig. 3 and agrees with 
our simulation result for a = 1. 

The above conclusion about the strong inaccuracy of 
the first Sonine approximation for a < 0.7 is further 
confirmed by Fig. 4, which shows the a-dependence of 
the three first Sonine coefficients. We observe that b'l is 
much less sensitive to dissipation than 63 ^nd ^3- For 
not too large inelasticity (a > 0.7), the second and third 
Sonine coefficients are much smaller than b'l and, conse- 
quently, one may expect that the series (2.44) truncated 
after fc = 1 is sufficiently close to the true distribution A! , 
at least in the region of thermal velocities relevant for the 
evaluation of the heat flux. On the other hand, as dis- 
sipation increases beyond a w 0.7, the magnitude of the 
coefficients b'2 and feg grow rapidly, becoming comparable 
to b'l. This implies that contributions to A beyond the 
first Sonine term cannot be neglected. In principle, the 
discrepancies of the first Sonine approximation could be 
remedied by considering the second Sonine approxima- 
tion, i.e., by truncating the expansion (2.44) after fc = 2. 
However, according to Fig. 4, there does not exist a range 




FIG. 3: (Color online) Plot of the modified thermal conduc- 
tivity coefficient n' versus the coefficient of normal restitution 
a. The circles refer to the simulation data obtained by the 
method described in this paper, while the triangles refer to the 
simulation data obtained in Ref. [11] from the GK relations. 
We have checked that the error bars of our simulation data 
are smaller than the size of the symbols. The solid line is the 
prediction given by the first Sonine approximation. The ar- 
row indicates the value k/kq = 1.025218 corresponding to the 
elastic limit when higher order terms in the Sonine expansion 
are taken into account. 
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FIG. 4: (Color online) Plot of the Sonine coefficients 6i (cir- 
cles), 62 (squares), and 63 (triangles) versus the coefficient of 
normal restitution a. The error bars are smaller than the size 
of the symbols. The lines are guides to the eye. 



of values of a where the magnitude of 63 is much smaller 
than that of 62, so that the second Sonine approximation 
possibly would not be sufficient. 

It is worthwhile remarking that the characteristic value 
of a approximately below which the first Sonine approxi- 
mation begins to deviate significantly from the simulation 
data is also the value below which the fourth cumulant 
a2 of the HCS starts to grow [32, 33, 38, 39]. This means 
that there seems to be a close relationship between the 
deviation of the HCS distribution /*^"^ from its Gaussian 
form and the deviation of the NS distribution /^^^ from 
its first Sonine approximation. Given that the weight 
function in the Sonine expansion is the Gaussian, it seems 
natural to conjecture that a better estimate of /'^' can 
be obtained by replacing the Gaussian weight function 
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by the HCS distribution /(°) [18, 28]. 

One of the advantages of the simulation method de- 
vised in this paper is that it allows to measure not only 
the transport coefficient k' or higher velocity moments 
(such as those associated with 63 a-nd 63), but also the 
NS velocity distribution function itself. As said in Sec. 
Ill, all the information contained in the full distribution 
j(i) (V) is present in the marginal distribution function 
g'^^^iYx) defined by Eq. (3.13). Thus, in what follows, we 
restrict ourselves to this marginal distribution function, 
which can be computed more efficiently than /'-^^(V) in 
the simulations. In order to analyze g^^'{Vx), it is conve- 
nient to write it in the form 



i(i) 



{V.) 



-1^-1/2,-c 



■■Cxip{cl)e\ 



(4.3) 



where ^{c^) is a dimensionless isotropic function, inde- 
pendent of e* . Its Sonine expansion is 



^{cD 



fe=i 



KLi'^'' 



{ell 



(4.4) 



where, as proven in the Appendix, the coefficients 6'j, are 
the same as in the Sonine expansion (2.44). For the sake 
of comparison, it is convenient to introduce the truncated 
series 



fc=i 



(4.5) 



If the first Sonine approximation is reliable, this means 
that (/'(c^) — '/'i(c^) = b'i{^ — c^) in the region of thermal 
velocities (say c^ < 10). In addition, one would expect 
that an even better approximation would be obtained 
with ((52 (c^) and </?3(c^)- 

The function (p(c^) is plotted in Figs. 5 (for a = 1 and 
a = 0.9) and 6 (for a — 0.7, a — 0.5, and a — 0.3). 
The first three truncated polynomials ipp{c^), p = 1, 2, 3, 
obtained by using the simulation values of b' are also 
plotted. In Fig. 5 we observe that the first Sonine polyno- 
mial fi{c1) captures reasonably well the behavior of the 
true distribution </?(c^), although it underestimates the 
latter for c^ > 4. The addition of the second and third 
Sonine polynomials significantly improves the agreement 
with the true distribution. In fact, cps is practically in- 
distinguishable from ip. As the inelasticity increases, so 
does the deviation of f from the first Sonine polynomial, 
as shown in Fig. 6. In contrast to the situation in Fig. 
5, ifi overestimates the function ip for large velocities. 
Moreover, the second and third Sonine truncated expan- 
sions do not clearly improve the agreement, especially for 
a = 0.5 and a — 0.3. 

The comparison carried out in Figs. 5 and 6 confirms 
that, ii a < 0.7, the velocity distribution function is not 
sufficiently well described by the first term in the Sonine 
expansion and, consequently, the value of b'l (and hence 
of k') estimated from the first Sonine approximation is 
not accurate. It is also interesting to note that when ipi 




FIG. 5: (Color online) Plot of the distribution function v'(c^) 
(solid lines) for a = 1 (top panel) and a = 0.9 (bottom panel). 
The dotted, dotted-dashed, and dashed lines correspond to 
the truncated Sonine expansions (pp{c^) with p = 1, 2, and 3, 
respectively. 



underestimates (overestimates) the function ip for large 
velocities, k' tends to be underestimated (overestimated) 
by the first Sonine approximation. In view of the panels 
corresponding to a = 0.5 and a = 0.3 in Fig. 6, it is 
highly questionable that the second or third Sonine ap- 
proximations could provide accurate estimates for k' . An 
alternative route has been recently proposed [28]. 

Finally, we address the question on the velocity range 
relevant to the evaluation of k' cx b'^. According to Figs. 
5 and 6, if that range were, for instance, < c^ < 6, 
then the first Sonine estimate would be reliable even for 
a = 0.3. Since this is not the case, it is obvious that the 
population of particles with c^ > 6 must have a signifi- 
cant influence on the heat flux for large inelasticity. To 
analyze this point with more detail, let us introduce the 
function 



Piicl) 



8 



37rV2 



du^ 



■^p{ul). (4.6) 



This function represents the contribution to b'^ coming 
from particles whose x-component of the (reduced) veloc- 
ity is smaller than \cx\, so that \iTa^2^ac /5i(c^) = b'^. The 
ratio Pi{Cx)/b'i is plotted in Fig. 7 for the same values of 
a as in Figs. 5 and 6. We see that the larger the inelastic- 
ity, the wider the range of velocities needed to faithfully 
evaluate b[. For instance, particles with < c^ < 6 con- 
tribute to 87% of the heat flux in the elastic case, while 
that percentage reduces to 66% in the case of a = 0.3. 
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FIG. 6: (Color online) Plot of the distribution function (p{c^) 
(solid lines) for a = 0.7 (top panel), a — 0.5 (middle panel), 
and a = 0.3 (bottom panel). The dotted, dotted-dashed, and 
dashed lines correspond to the truncated Sonine expansions 
Vp(c^) with p = 1, 2, and 3, respectively. 



V. CONCLUSIONS 

The CE solution of the inelastic Boltzmann equation 
provides expressions for the NS transport coefficients in 
terms of the solutions of linear integral equations [13] . Al- 
ternatively, these expressions can be proven to be equiva- 
lent to GK relations [5]. In either case, in order to obtain 
the explicit dependence of the transport coefficients on 
the coefficient of restitution, one needs to make use of 
certain approximations. As in the case of elastic colli- 
sions, the standard approach is to expand the unknown 
NS velocity distribution function in Sonine polynomials 
and truncate it at a given order. Of course, the practi- 
cal difficulties increase considerably with the number of 
retained polynomials, so the first Sonine approximation 
is usually chosen. To check the reliability of the first 
Sonine approximation one needs to resort to compari- 
son with DSMC computer simulations of the Boltzmann 
equation. In principle, there are two strategies to mea- 
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FIG. 7; (Color online) Plot of the ratio /3i(c^)/&i for a = 1, 



sure the transport coefficients via computer simulations 
in homogeneous states. One consists of measuring the ap- 
propriate time correlation functions in the HCS and then 
carry out an integration over time by applying the GK 
relations [5]. In the other possible strategy, one weakly 
disturbs the granular gas from the HCS by the action of 
a conveniently chosen homogeneous, anisotropic external 
force that produces the same effect as a hydrodynamic 
gradient; in this way, the simulations provide in the lin- 
ear response regime the transport coefficients as well as 
the NS velocity distribution functions. The first method 
has been used to get the shear viscosity r\ and the two 
transport coefficients (k and /i) associated with the heat 
flux [10, 11], while so far the second method has only 
been applied to r\ [21]. 

The first Sonine approximation for r\ is seen to com- 
pare quite well with DSMC simulations, even for strong 
dissipation [10, 11, 21]. However, the corresponding ap- 
proximations for K and /i appreciably differ from GK 
simulations for high inelasticity [a < 0.7). One could 
argue that these discrepancies are a reflection of the ve- 
locity correlations appearing in the HCS, even in the low- 
density limit. In that scenario, the correlation functions 
computed in DSMC simulations would incorporate effects 
not accounted for by the Boltzmann equation. 

The primary goal of this paper has been to investigate 
whether the discrepancies observed between the first So- 
nine estimates and the GK data are due to effects beyond 
the Boltzmann framework or are simply due to the limita- 
tions of the first Sonine approach. To that end, we have 
followed the second strategy mentioned above, namely, 
the one based on the action of a homogeneous external 
force. First, we have derived the explicit forms of the re- 
spective velocity-dependent external forces which yield, 
to ffist order, the NS velocity distribution functions asso- 
ciated with the standard thermal conductivity k and the 



modified thermal conductivity k! 



n^./2T. In the 



case of K, the external force must be complemented by a 
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stochastic term describing a donation-annihilation pro- 
cess. Since the presence of this latter term complicates 
the simulation method, in this paper we have focused 
on the homogeneous Boltzmann equation corresponding 
to the coefficient k' . Our DSMC results show an ex- 
cellent consistency with those obtained in Ref. [11] by 
the GK formalism. Since our method is entirely tied to 
the Boltzmann equation, we conclude that the deviations 
of the simulation data from the first Sonine approxima- 
tion are mainly due to the inaccuracy of the latter. This 
conclusion has been further supported by an analysis of 
the first three Sonine coefficients and of the NS velocity 
distribution function itself. While the second and third 
Sonine coefhcients are practically negligible for a > 0.7, 
they rapidly increase in magnitude as the inelasticity in- 
creases, becoming comparable to the first Sonine coeffi- 
cient. With respect to the distribution function, we have 
observed that, for strong dissipation, it is not well cap- 
tured by the first Sonine polynomial in the whole velocity 
region relevant to the computation of the transport coef- 
ficient k' . 

The simulation results reported in this paper indicate 
that the second or third Sonine approximations are not 
expected to improve significantly the quality of the first 
Sonine estimate, especially considering the technical dif- 
ficulties associated with the method. An alternative av- 
enue under the form of a modified first Sonine approxi- 
mation is explored elsewhere [28]. 
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APPENDIX: PROOF OF THE EQUIVALENCE 
BETWEEN EQS. (3.10) AND (3.12) 

The full NS velocity distribution function (3.8) can be 
written as 

/(i)(V) = nV''7r-'^/2e-^'c,<&(c2)e*, (A.l) 

where $(c^) is a dimensionless isotropic function. This 
function can be expanded in Laguerre (or Sonine) poly- 
nomials as 



<i>(c^)=5:5',Lr)(c^), 



(A.2) 



fc=l 



in agreement with Eq. (2.44). The orthogonality condi- 
tion (2.12) yields 



b'. 



^il±^l^n-'/^ fdccle-^'Lf'\c'Mc'). 



r{k + i + d/2) 



This expression is equivalent to Eq. (3.10). 

Let us now consider the marginal distribution g^^'{Vx) 
defined by Eq. (3.13). Similarly to Eq. (A.l), g'-^>{Vx) 
can be expressed as Eq. (4.3). Thus, the function (^(c^) 
can be obtained from <l>(c^) as 



^(c2) = ^-(-^-D/^ /dcxe-='^$(c2) 



(A.4) 



Inserting the expansion (A.2) we get 

oo 

(^(c,)=^6',Ffc(c2), 
fc=i 

where we have called 



(A.5) 



(A.6) 



Now we make use of the mathematical property [31] 



and take p — 1/2, q = (d — 3)/2, x = c^, and y = c\. 
Thus, 



(A.S 
The integral can be computed as 



.-(-i)/^/dc,e-14T^(4) 



1 



r{{d-i)/2)J, 
— Skj, 



dyy(^-^y^L['f\y) 



(A.9) 



where in the last step we have made use of the orthog- 
onality relation of the generalized Laguerre polynomials. 
Therefore, i^fc(c^) = ^^ (c^) and so Eq. (A.S) becomes 



^icl)^Y.b',L^^^'^{cl). (A.IO) 



From here one can obtain an alternative expression for 
fficients b'j, as 

2r(3/2)fc! _i. 



the coefficients b'j, as 



(A.3) 



(A.l 
This is equivalent to Eq. (3.12). 
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